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ABSTRACT 

We  use  the  commercial  computational  fluid  dynamics  code  Fluent  to  simulate  the  flow 
around  a  sphere  in  several  different  flow  regimes;  steady-state  laminar  flow  at  a  Reynolds 
number  (Re)  of  100,  time-dependent  laminar  flow  at  Re  =  300,  and  turbulent  flow  at  Re  =  104 
and  Re  =  106.  These  simulations  provide  a  test  of  the  ability  of  the  code  to  accurately 
reproduce  typical  flow  structures  observed  in  generic  bluff  body  flows,  such  as  those 
experienced  by  submarines  and  Unmanned  Underwater  Vehicles  (UU  Vs).  The  simulations  are 
compared  both  with  experimental  results  and  computations  from  other  computer  codes  and  it 
is  found  that  Fluent  is  able  to  accurately  simulate  the  fluid  behaviour  in  each  of  the  above 
flow  regimes. 
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Simulation  of  Flow  Past  a  Sphere 
usi  ng  the  FI  uent  Code 


Executi  ve  Su  mmary 

The  use  of  Computational  Fluid  Dynamics  (CFD)  codes  to  simulate  the  flow  around 
geometrically  complicated  shapes  such  as  aeroplanes,  cars  and  ships  has  become 
standard  engineering  practice  in  the  last  few  years.  The  computer  code  Fluent  is  a 
commercial  CFD  code  which  is  used  routinely  by  members  of  the  Flydrodynamics 
Group  within  DSTO's  Maritime  Platforms  Division  (MPD)  to  simulate  flows  relevant 
to  underwater  vehicle  design  and  surface  ship  wakes.  In  this  report  we  test  the  ability 
of  the  code  to  accurately  simulate  the  flow  around  a  sphere  over  a  large  range  of 
Reynolds  numbers  (Re).  This  is  a  stringent  test  of  the  code  as  it  involves  a  number  of 
very  different  flow  regimes,  varying  from  laminar  steady-state  flow  near  Re  =  100  to 
time-dependent  turbulent  flow  at  Re  =  106. 

These  simulations  are  also  intended  to  provide  a  benchmark  comparison  for  the  finite 
difference  CFD  code  Vortel,  currently  under  development  at  Naval  Underwater 
Warfare  Center,  Newport,  Rhode  Island  in  collaboration  with  researchers  at  MPD. 
Vortel  is  based  on  the  Lagrangian  Vorticity  method  and  offers  advantages  over 
commercial  CFD  codes  for  the  simulation  of  flow  around  multiple  bodies  in  relative 
motion.  It  has  already  been  used  to  simulate  unmanned  undersea  vehicle  docking 
manoeuvres,  as  well  as  unsteady  bow  thruster  hydrodynamics  and  the  unsteady 
separated  flow  fields  past  an  oscillating  airfoil. 

Four  different  flow  regimes  are  studied  in  detail:  steady-state  laminar  flow  at  Re  =  100, 
time-dependent  laminar  flow  at  Re  =  300,  turbulent  flow  with  laminar  boundary  layers 
at  Re  =  104  and  turbulent  flow  with  turbulent  boundary  layers  at  Re  =  106.  The 
simulated  flows  were  found  to  be  in  excellent  agreement  with  both  experimental 
results,  where  available,  and  with  the  results  obtained  by  other  authors  using  different 
simulation  codes.  The  simulations  described  in  this  report  illustrate  the  capability  of 
the  Fluent  code  to  accurately  reproduce  typical  flow  structures  observed  on  this 
generic  bluff  body  flow  for  both  time  independent/  time-dependent  and 
laminar/ turbulent  flow  regimes. 
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1.  Introduction 

The  use  of  computational  fluid  dynamics  codes  to  simulate  the  flow  around  geometrically 
complicated  shapes  such  as  aeroplanes,  cars  and  ships  has  become  standard  engineering 
practice  in  the  last  few  years.  A  number  of  commercially  available  codes  can  be  used  to 
perform  these  studies,  including  the  finite  volume  codes  Fluent  [11]  and  CFX  [3],  both  of 
which  are  used  routinely  by  the  FI  ydrodynamics  Group  within  DSTO's  Mari  time  Platforms 
Division  (M  PD)  to  perform  computational  fluid  dynamics  (CFD)  studies  on  underwater  flow 
problems  relevant  to  underwater  vehicle  design.  The  strength  of  the  commercial  codes 
generally  lies  in  their  excellent  pre-processor  software,  which  allows  robust  meshes  to  be 
constructed  around  very  complicated  shapes,  and  the  number  of  quite  advanced  turbulence 
models  contained  within  them.  Whilst  these  codes  are  now  also  developing  the  ability  to 
handl  e  probl  ems  with  movi  ng  and  deformi  ng  meshes,  these  features  aretypi  cal  ly  structured 
around  particular  geometries  of  commercial  importance,  such  as  piston  motion  in  cylinder 
bl  ocks,  fan  bl  ades  i  n  cyd  ones  and  stores  separati  on  from  ai  rcraft.  The  abi  I  ity  to  si  mul  atefl  ow 
around  an  arbitrary  number  of  bodies  in  relative  motion,  which  is  of  interest  to  the  MPD 
FI  ydrodynamics  Group,  is  not  intrinsically  suited  to  these  codes. 

Onecomputercodewhich  shows  consi  derablepromi  sein  this  area  is  Vortel,  a  finite  difference 
codecurrently  under  development  at  NUWC,  Newport,  Rhodelsland  in  collaboration  with 
researchers  at  M  PD  [15].  Vortel  isbased  ontheLagrangian  Vorti  city  method  [35]  andwill  be 
used  to  simulate  flow  around  multi  plebodies  in  relative  motion.  It  has  already  been  used  to 
si  mul  ate  unmanned  undersea  vehicle  docking  manoeuvres  [17,18]  aswell  asunsteady  bow 
thruster  hydrodynami  cs  [  16]  and  the  unsteady  separated  fl  ow  fi  el  ds  past  an  osci  1 1  ati  ng  ai  rfoi  I 
[19],  As  part  of  the  development  process  Vortel  i  s  bei  ng  tested  agai  nst  a  number  of  benchmark 
problems,  one  of  which  is  the  simulation  of  the  flow  past  a  sphere  over  a  large  range  of 
Reynol  ds  numbers  (Re).  Thi  s  is  a  stri  ngent  test  of  the  code  as  it  i  nvolves  a  number  of  different 
fl  ow  regi  mes,  from  I  ami  nar  steady- state f  I  ow  near  Re = 100  to  ti  me-dependent  turbul  ent  fl  ow 
at  Re  =  106.  To  provide  verification,  the  Fluent  code  has  been  used  to  simulatetheflow  pasta 
sphere  in  several  flow  regimes;  steady-state  I  ami  nar  flow  at  Re =100,  time-dependent  I  ami  nar 
flow  at  Re  =  300,  and  turbulent  flow  at  Re  =  l<f  and  Re  =  106. 

As  Fluent  is  used  to  perform  the  majority  of  the  flow  simulation  studies  performed  by  the 
Fl  ydrodynamics  Group,  the  computations  reported  here  also  provide  a  timely  test  of  the 
ability  of  thecodeto  accurately  reproduce  typical  flow  structures  observed  in  generic  bluff 
body  flows.  Thesimulationsarecompared  both  with  experimental  results  and  computations 
f  rom  oth er  comp  u ter  cod  es  an d  i  t  i  s  f ou  n d  th at  F I  u  ent  i  s  abl  e  to  accu  ratel  y  si  mu  I  ate  th e  f  I  u  i  d 
behaviour  in  each  of  the  above  flow  regimes.  At  Re  =  300  the  laminar  solver  in  Fluent 
accurately  calculates  the  regular  vortex  shed  ding  observed  both  experimentally  and  in  many 
numerical  simulations,  while  the  new  Large  Eddy  Simulation  (LES)  turbulence  model 
avai  I  abl  e  i  n  th  e  I  atest  versi  on  of  F I  u  ent  ( F I  u  ent  6. 3)  i  s  sh  ow  n  to  be  abl  e  to  si  mu  I  ate  th  e  d  rop  i  n 
drag  coefficient  which  occurs  near  Re  =  3.8  x  1&  due  to  the  transition  from  laminar  to 
turbulent  boundary  layer  separation.  Comparison  with  the  simulation  results  obtained  from 
the  Vortel  code  are  described  elsewhere  [15]. 
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2.  Flow  regimes  as  a  function  of  Reynolds  number 

The  nature  of  the  flow  around  a  sphere  changes  dramatically  as  the  Reynolds  number  of  the 
flow  increases.  In  general,  the  higher  the  Reynolds  number,  the  morecomplextheflow.  I  n  this 
secti  on  w  e  d  escri  be  the  natu  re  of  the  f  I  ow  for  d  i  fferent  Reynoldsnu  mbers  i  n  a  range  betw  een 
Re  =  20and  Re  =  1.0xl06.  Itshould  be  noted  that  the  Reynolds  number  range  for  each  of  the 
different  fl  ow  regi  mes  shows  some  vari  ati  on  dependi  ng  on  the  i  ndi  vi  dual  researchers,  so  that 
the  Reynolds  numbers  given  here  delineating  one  flow  regime  from  another  should  be 
considered  to  havea  small  amount  of  uncertainty  attached  to  them. 

2.1  Steady  axisymmetric  regime:  20  <  Re  <  210 

Experimental  investigations  of  the  steady  wake  behind  a  sphere  at  low  Reynolds  numbers 
have  been  performed  byTaneda[48]  and  Nakamura  [38],  Taneda  found  that  for  Reynolds 
numbers  less  than  24  the  flow  around  the  sphere  is  perfectly  laminar,  no  flow  separation 
occurs,  and  theflow  onthedownstream  sideofthesphereis  identical  to  that  on  the  upstream 
side. 

Above  a  Reynolds  number  of  25  the  flow  separates  from  the  sphere  close  to  the  rear 
stagnati  on  poi  nt  and  forms  a  cl  osed  red  rcu  I  ati  ng  w  ake  i  n  the  shape  of  an  axi  sy  mmetri  c  vortex 
ri  ng.  A  s  the  Rey  nol  ds  nu  mber  i  ncreases  both  the  separati  on  angl  e  and  the  I  ength  of  the  wake 
grow.  Taneda  showed  that  the  size  of  the  vortex  ring  was  proportional  to  the  logarithm  of  the 
Reynolds  number  and  that  theflow  was  stationary  up  to  a  Reynolds  number  of  approximately 
130,  after  w  hi  ch  osci  1 1  ati  ons  began  to  appear  at  the  rear  of  the  vortex  ri  ng.  N  akamu  ra  however 
found  that  the  vortex  ring  was  stable  and  axisymmetric  up  to  a  Reynolds  number  of 
approxi  mately  190,  and  that  a  stabl  evortex  ri  ng  formed  at  a  Reynolds  number  as  I  ow  as  Re  = 
7. 

M  agarvey  and  Bishop  [34]  and  Wu  and  Faeth  [54]  experimentally  investigated  theflow  pasta 
sphere  over  a  larger  range  of  Reynolds  numbers  and  found  similar  behaviour  to  that  found  by 
Nakamura  [38],  Wu  and  Faeth  found  that  theflow  was  axisymmetric  and  stable  up  to  Re  = 
200,  while  M  agarvey  and  Bishop  found  the  same  behaviour  occurring  up  to  Re  =  210.  These 
observations  are  in  good  agreement  with  the  calculations  of  Natarajan  and  Acrivos[39],  who 
investigated  the  I  inear  stability  of  the  steady  axi  symmetric  flow  past  a  sphere  and  found  that 
theflow  undergoes  a  regul  ar  bifurcation  at  a  Reynolds  number  of  about  210 and  results  i  n  the 
development  of  a  non-axi symmetric  wake. 

2.2  Steady  planar-symmetric  regime:  210 <  Re <  270 

Above  this  transition  point  near  Re  =  210  M  agarvey  and  Bishop  [34]  and  Wu  and  Faeth  [54] 
noted  that  theflow  remained  attached  and  stable  but  was  no  longer  axi  symmetric.  Thenature 
of  theflow  in  this  regime  consists  of  two  streamwise  vortical  tails  of  equal  strength  and 
opposite  sign.  M  agarvey  and  Bishop  [34]  referred  to  this  structureas  a  doublethreaded  wake 
and  were  able  to  observe  and  photograph  the  phenomenon  by  following  the  drop  of  an 
immisciblefluid  in  a  number  of  different  liquid-liquid  systems. 


2 


DSTO-TR-2232 


Althoughtheflow  no  longer  possesses  axial  symmetry  theflow  still  exhibits  planar  symmetry 
i  n  the  pi  ane  contai  ni  ng  the  two  vorti  cal  tai  I  s.  The  ori  entati  on  of  thi  s  symmetry  pi  ane  i  s  ti  I  ted 
off  the  streamwise  axis  and  has  been  explained  byjohnson  and  Patel  [23]  as  being  dueto  an 
azimuthal  instability  of  the  low  pressurecoreofthewakevortex.Theactual  orientation  of  the 
plane  is  arbitrary  and  will  be  determined  by  random  external  influences,  such  as 
perturbations  dueto  model  supports  in  an  experimental  situation,  or  truncation  errors  or  grid 
asymmetries  in  numerical  simulations. 

The  I  oss  of  axi  al  symmetry  results  i  n  the  appearance  of  a  I  ateral  (I  ift)  force  whi  ch  i  s  di  rected  i  n 
the  pi  ane  of  flow  symmetry.  Thepresenceofthisforcewas  observed  by  M  agarvey  and  Bishop 
[34]  by  noting  that  the  liquid  drops  in  their  experiments  deviated  from  a  vertical  line  of  fall. 
More  recently  the  presence  of  a  lateral  force  has  been  verified  by  several  numerical 
simulations.  Gushchin  et  al.  [14]  used  an  explicit  second  order  finite  difference  scheme  to 
perform  d  i  red  numeri  cal  si  mu  I  ati  ons  of  thef  I  ow  arou  nd  a  sphere  and  fou  nd  zero  I  ift  force  at 
Re  =  210  but  a  finite  value  at  Re  =  210.5.  Johnson  and  Patel  [23]  found  that  the  simulated  lift 
force  jumped  by  three  orders  of  magnitude  between  Re  =  211  and  Re  =  212. 

Tomboulides  and  Orszag  [50]  used  dired  numerical  simulation  based  on  spedral-type 
methodsto  simulate  theflow  between  Re=25and  Re  =1000.  Their  simulations  showed  that 
theflow  past  a  sphere  is  axi  symmetric  up  to  a  Reynolds  number  of  approximately  212,  and 
that  beyond  this  Reynolds  number  theflow  undergoes  a  transition  to  three-dimensionality 
th  rou  gh  a  regu  I  ar  bi  f  u  rcati  on .  The  charaderi  sti  cs  of  the  resu  I  ti  ng  steady  f  I  ow  f  i  el  d ,  w  hi  ch  they 
found  to  bestableup  to  a  Reynolds  number  of  approximately  270,  corresponded  to  that  of  the 
double-thread  wake  reported  by  M  agarvey  and  Bishop  [34].  Thesimulati  ons  showed  that  the 
double  thread  wake  consisted  of  two  opposite-sign  streamwise  vortices  which  extend  to 
infinity  and  appear  in  the  experiment  as  two  dye  threads  emanating  from  the  end  of  the 
recirculation  region. 

2.3  U nsteady  planar-symmetric  regime:  270 <  Re <  400 

As  the  Reynolds  number  increases  within  this  range  a  transition  from  the  steady  planar- 
symmetri  c  wake  to  a  ti  me-dependent  planar-symmetric  wake  occurs.  The  unstead  i  ness  fi  rst 
appears  as  a  wavi ness  i  n  thedoubl e threaded  wake.  Thi sform  of  i nstabi lity  in  the wakeexi sts 
only  over  a  very  limited  large  of  Reynolds  numbers,  between  Re  =  270  and  Re  =  290.  Above 
Re = 290  M  agarvey  and  Bishop  [34]  found  that  the  wake  became  unsteady  and  consisted  ofa 
series  of  interconnected  vortex  loops,  while  Wu  and  Faeth  [54]  found  that  theflow  became 
unsteady  and  exhibited  vortex  shedding  at  a  Reynolds  number  close  to  280.  Sakamoto  and 
Haniu  [43]  conducted  an  extensive  series  of  experiments  on  the  flow  past  spheres  in  the 
Reynolds  number  range  from  300  to  40,000  using  hot-wire  techniques  in  a  low  speed  wind 
tunnel  and  fl  ow-vi  sual  i  sati  on  techni  ques  i  n  a  water  channel .  They  found  that  hai  rpi  n  shaped 
vortices  began  to  be  shed  periodically  at  Re  =  300  with  regularity  in  strength  and  frequency 
and  that  the  planar  symmetry  of  theflow  was  maintained. 

These  experimental  observations  [34, 43, 54]  arein  broad  agreementwith  a  number  of  recent 
numerical  simulations.  Johnson  and  Patel  [23]  performed  numerical  simulations  at  Re  =270 
and  found  a  steady-state  solution,  while  at  Re  =  280,  the  next  highest  Reynolds  number 
considered  in  their  calculations,  the  solution  clearly  became  periodic.  The  results  of  their 
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simulation  at  Re  =  300  showed  a  highly  organised  periodic  flow  dominated  by  vortex 
shedding,  with  a  Strouhal  number  of  0.137  and  a  mean  valuefor  the  drag  coefficient  Cd  of 
0.656. 

Ploumhans,  et  al.  [40]  used  three-dimensional  vortex  methods  to  si  mu  I  ate  the  flow  past  a 
sphere  at  Reynolds  numbers  of  300,  500  and  1000.  At  Re  =  300  they  found  that  the  flow 
attained  a  time  periodic  regimewith  a  Strouhal  number  of  0.135  and  amean  valuefor  Cd  of 
0.683.  A  finite  value  for  the  lift  coefficient  was  obtained  and  the  flow  exhibited  planar 
symmetry  i  n  a  plane  passi  ng  through  the  wake  centre  I  i  ne. 

The  direct  numerical  simulations  of  Tomboulides  and  Orszag  [52]  used  a  mixed  spectral 
element/  Fourier  spectral  method  and  found  that  the  transition  between  three-dimensional 
steady  flow  with  planar  symmetry  to  single  frequency  periodic  flow  with  vortex  shedding 
occurred  in  the  range  between  Re  =  270  and  Re  =  285.  They  also  performed  a  simulation  at 
Re  =  300  and  found  a  Strouhal  number  of  0.136  and  a  mean  drag  coefficient  of  0.671.  They 
noted  that  the  basic  wake  structure  consisted  of  a  succession  of  interconnected  vortex  loops 
which  maintained  the  planar  symmetry  observed  at  lower  Reynolds  numbers. 

Mittal  [37]  performed  numerical  simulations  of  flow  past  prolate  spheroids  using  a  Fourier- 
Chebyshev  spectral  collocation  method  for  Reynolds  numbers  in  the  range  from  Re  =  50to 
Re  =  500.  At  Re  =  350  he  found  that  the  wake  was  characterised  by  the  appearance  of 
interconnected  vortex  loops  which  were  shed  from  the  sphere  at  a  well  defined  frequency 
with  a  Strouhal  number  of  0.140.  Planar  symmetry  could  also  be  observed  in  the  unsteady 
wake  in  a  plane  passing  through  the  wake  centre  line. 

2.4  U nsteady  asymmetric  regime:  400 <  Re <  1000 

I  n  thi  s  range  the  pi  anar  symmetry  of  thefl  ow  i  s  eventual  ly  lost  because  the  azi  muthal  I  ocati  on 
at  w  hi  ch  the  vortex  I  oops  are  formed  changes  from  cycl  e  to  cyd  e  i  n  an  i  rregu  I  ar  fashi  on .  The 
exact  value  of  the  Reynolds  number  at  which  this  begins  to  happen  is  difficult  to  specify 
precisely.  In  the  experiments  of  Sakamoto  and  Flaniu  [43]  it  was  noted  that  the  shedding  of 
the vorti  ces  started  to  become  i  rregul ar  near  Re  =420 and  was  compl  etely  random  at  Re  =480. 
Mittal  [36]  conducted  a  series  of  direct  numerical  simulations  in  the  range  350  <  Re<425to 
more  accurately  locate  the  upper  extent  of  the  planar  symmetric  range  and  found  that  the 
wake  I  oses  pi  anar  symmetry  betw  een  Re  =  350  and  Re  =  375.  The  cycl  e-to-cycl  e  vari  ati  ons  i  n 
thevortex  formation  angleattheseReynoldsnumberswherequitesmall  however  and  did  not 
becomeappreciableuntil  Re=425.  Mittal  [36]  thus  explained  this  apparent  discrepancy  with 
the  results  of  Sakamoto  and  Flaniu  as  being  due  to  the  inability  of  experimental  flow 
visualisation  techniques  to  detect  small  scale  variations  in  the  azimuthal  angle  of  vortex 
formation.  The  numerical  simulations  of  Dallmann  et  al.  [9]  however  indicate  that  planar 
symmetry  still  exists  in  the  flow  at  R  =400. 

A  t  Re  = 500 the  si  mu  I  ati  ons  of  T ombou  I  i  d  es  and  O  rszag  [  50]  show  ed  a  w  ake  structu  re  si  mi  I  ar 
to  the  one  observed  at  Re =300,  but  thevortex  I  oops  were  shed  from  thespherewith  different 
orientations,  hence  the  planar  symmetry  observed  at  Re  =  300  was  not  preserved.  Time 
spectra  of  the  velocity  at  several  locations  showed  that  the  flow  had  changed  from  a  one- 
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frequency  flow  at  Re =300 to  an  al  most  chaoti  c  system  at  Re =500,  although  a  domi  nant  peak 
at  a  Strouhal  number  of  0.167  could  still  be  observed. 

A  further  transition  in  the  nature  of  the  flow  occurs  near  Re  =  800.  Above  this  Reynolds 
number  the  measurements  of  Kim  and  Durbin  [26]  showed  that  two  dominant  modes  of 
unsteadiness  exist  in  the  wake.  These  are  associated  with  a  small  scale  instability  in  the 
separating  shear  layer  (the  Kelvin-Helmhotz  instability)  and  a  large  scale  instability  in  the 
w  ake  (the  vortex  shedd  i  ng).  The  hi  gher  frequency  component  was  onl  y  detected  i  n  the  regi  on 
immediately  downstream  of  the  sphere,  leading  Kim  and  Durbin  to  the  conclusion  that  this 
modewas  dueto  a  Kelvin-Helmhotz  instability  in  the  separating  shear  layer.  The  Strouhal 
number  associated  with  vortex  shed  ding  was  found  to  be  practically  independent  of  Reynolds 
number  for  Re  >  800  and  had  a  value  of  approximately  0.2,  while  the  Strouhal  number 
associated  with  the  Kelvin- Helmholtz  instability  increased  with  increasing  Reynolds  number. 
Similar  results  have  also  been  found  in  the  experiments  of  Sakamoto  and  Haniu  [43]  and 
C  homaz  et  al .  [4] .  The  presence  of  a  Kel  vi  n-  H  el  mhol  tz  i  nstabi  I  ity  i  n  the  separati  ng  shear  I  ayer 
above  Re  =800  has  also  been  seen  in  the  numerical  simulations  of  Tomboul  ides  and  Orszag 
[50],  who  used  a  Spectral  Element/  Fourier  Spectral  method. 

2.5  Turbulent  wake  regime:  Re  >  1000 

T  aneda  [49]  made  vi  sual  observati  ons  of  thefl  ow  past  a  sphere  i  n  a  wi  nd  tunnel  i  n  the  range 
104<  Re<  106  using  a  combi  nation  ofthesurfaceoil  flow  method,  the  smoke  flow  visualisation 
method  andthetuft-grid  method.  For  Reynolds  numbers  in  therangel.Ox  104<  Re<  3.8x10^ 
he  showed  that  the  wake  undergoes  a  progressive  wave  motion  in  a  plane  containing  the 
streamwise  axis.  The  Strouhal  number  and  drag  coefficient  remain  approximately  constant  at 
0. 2  and  0. 5  respecti  vel  y ,  but  the  pi  ane  contai  n  i  ng  the  osci  1 1  ati  ons  rotates  si  ow  I  y  and  i  rregu  I  arl  y 
around  the  axis.  These  observations  are  in  agreement  with  those  of  Achenbach  [1,2],  who 
found  that  strong  periodic  fluctuations  existed  in  thewake  right  up  to  the  critical  Reynolds 
number  of  3  x  1&,  and  that  the  vortex  separation  point  rotated  around  the  sphere. 

Above  the  critical  Reynolds  number  around  Re=3.8x  1&  the  near-wake  recirculation  region 
shrinks  considerably,  the  drag  coefficient  decreases  sharply  to  a  value  of  Cd  «  0.08,  and 
peri  od i  c  vortex  shedd i  ng  i  s  no  I  onger  abl  eto  be  detected  experi  mental  I y .  The  experi  ments  of 
Taneda[49]  show  that  in  the  range  3.8  x  l(f<  Re<  1.0  x  106  the  sphere  wake  forms  a  pair  of 
streamwise  line  vortices  at  a  short  distance  from  the  streamwise  axis  and  the  vortex  pair 
rotates  slowly  and  randomly  about  that  axis. 

Numerical  si  mu  I  ati  ons  of  flow  for  Reynolds  numbers  greater  than  103  have  been  performed 
mainly  using  LargeEddy  Simulation  (LES)  or  Detached  Eddy  Simulation  (DES)  techniques, 
although  a  few  authors  have  attempted  to  use  either  Reynolds  Averaged  Navier  Stokes 
(RANS)  or  unsteady  Reynolds  Averaged  Navier  Stokes  (URANS)  methods,  with  limited 
success.  Drikakis[10]  computed  the  flow  past  a  sphere  at  Reynolds  numbers  of  1.62  x  lCf  and 
5xl06using  both  ak-e  model  and  the Baldwin-Lomax  model.  Comparison  with  experiment 
was  I  imited  to  plots  of  pressurecoeffi  dent  versus  polar  angle.  Atthelower  Reynolds  number 
the  k- e  model  gave  good  agreement  with  the  experimental  results  prior  to  the  separation 
poi  nt,  but  I  i  mited  agreement  i  n  the  aft  regi  on.  At  the  hi  gher  Reynol  ds  number  the  k-s  model 
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performed  poorly  for  most  values  of  the  polar  angle.  The  Baldw in-Lomax  model  was  very 
inaccurate  at  both  Reynolds  numbers. 

Constantinescu  and  Squires  [6]  simulated  the  flow  over  a  sphere  at  Re  =  lCf  using  URANS, 
LESand  DEStechniques.TheURA  NS  calculations  were  performed  using  thek-s,/c-co,  V2f  and 
Spalart-Allmaras  models.  Spatial  discretisation  was  based  on  generalised  curvilinear 
coordi  nates  and  high  order  fi  nite  difference  approxi  mati ons  and  ti  me  advancement  was  vi a 
thefractional  step  method.  With  theexception  of  the  two-  layer  k-s  model,  which  performed 
poorly  everywhere,  the  URANS  predictions  of  the  pressure  coefficient  and  thestreamwise 
drag  coefficient  were  in  reasonable  agreement  with  experiment.  The  URANS  simulations 
however  were  not  able  to  accurately  simulate  the  drag  coefficient  and  the  energy  content 
associ ated  with  the  mai  n  i  nstabi  I  ity  modes.  Both  the  L ES  and  D ES  computati ons  showed  the 
presence  of  a  second  high  frequency  band  in  the  power  spectra  and  the  calculated  Strouhal 
numbers  for  both  the  low  frequency  (vortex  shedding)  and  high  frequency  (Kelvin-H  elmholtz 
instability)  modes  were  in  good  agreement  with  the  experimental  results  of  Sakomoto  and 
Haniu  [43].  None  of  the  URANS  simulations  [6,10]  showed  the  presence  of  the  higher 
frequency  band. 

Constantinescu  and  Squires  [7]  simulated  theflow  over  a  sphere  using  DESon  a  structured 
hexahedral  mesh  having  450,000  nodes  for  Reynolds  numbers  of  1.0  x  1C f,  4.2  x  lCf  and 
1.14x  106.  In  thesubcritical  regime  (Re  =  1.0  x  1(f)  with  laminar  boundary  layer  separation 
D  ES  was  ableto  capture  the  largescale  vortex  shedding  as  well  astheformation  of  the  vortex 
tubes  in  the  separated  shear  layers  due  to  the  growth  of  Kelvin-H  elmholtz  instabilities  in 
theselayers.  Good  agreement  was  found  with  experimental  data  for  drag  coefficient,  location 
of  transition  on  the  sphere,  and  the  distribution  of  pressure  and  skin-friction  coefficients.  In 
thesupercritical  regime(Re=4.2x  lCfandl.Mx  106)  DESaccurately  predicted  theposition of 
boundary  layer  separation  and  the  distribution  of  the  mean  pressure  coefficient  over  the 
sphere  surface,  but  the  skin  friction  coefficient  was  less  accurately  predicted.  Improved 
agreement  w  as  obtai  ned  w  hen  the  model  was  acti  vated  near  the  regi  on  w  heretransi  ti  on  w  as 
observed  in  experiments. 

Kim  [29]  used  LES  on  unstructured  meshes  to  simulate  flow  around  a  sphere  for  Reynolds 
numbers  of  1.0  x  104and  1.14  x  106.  Two  different  sub-grid  scale  stress  models  were  used,  a 
dynamic  Smagorinsky  model  and  a  localised  dynamic  kinetic  energy  model,  but  there  was 
little  difference  in  the  results  predicted  by  the  two  models.  At  Re  =  1.0  x  104both  models 
accurately  predicted  the  magnitude  of  the  drag  coefficient,  the  vortex  shedding  Strouhal 
number  and  the  separation  angle.  As  Re  increased  further  the  magnitude  of  the  drag 
coefficient  was  accurately  predicted  but  the  locations  of  the  separation  point  and  the  onset  of 
the  transition  from  laminar  to  turbulent  flow  in  the  boundary  layer  were  much  less 
satisfactory.  A  summary  of  the  different  flow  regimes  is  contained  inTablel. 
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Table  1:  Summary  of  different  flow  regimes  as  a  function  of  Reynolds  number 


Re  range 

Nature  of  Flow 

Drag  coefficient 

Strouhal  number 

20<  Re<  210 

Steady  axi  symmetric  flow  with  closed 
red  rcu  1  ati  ng  wake  i  n  the  shape  of  a  vortex  ri  ng. 

2.7  — »  0.7 

Steady  state  flow 

210  <  Re  <270 

Flow  is  attached,  stable  and  consists  of  two 
streamwise  vortical  tails  of  equal  strength  and 
opposite  sign.  Flow  exhibits  planar  symmetry. 

Cd  continues  to  fall 
monoton  i  cal  ly 

Steady  state  flow 

270  <  Re  <300 

Transition  to  unsteady  wake  consisting  of  the 
peri  od  i  c  shedd  i  ng  of  a  seri  es  of  i  ntercon  nected 
vortex  loops.  Flow  still  exhibits  planar 
symmetry. 

Mean  value  of 

Cd®  0.66 -0.68 

St  *  0.14 

400  <  Re  <800 

Vortex  loops  continue  to  be  shed  but  the 
azimuthal  angle  at  which  they  are  formed 
begins  to  oscillate  in  an  irregular  fashion  and 
the  planar  symmetry  is  lost. 

Mean  valuefallsto 

0.5  at  Re  =  800 

St  slowly  increases 
to  0.2  at  Re  =  800 

800  <  Re  <  3.8  x  105 

Vortex  loops  continue  to  be  shed  but  in  a  more 
irregular  manner,  although  strong  periodic 
fluctuations  still  exist  in  the  wake. 

Cd  »  0.5  over  entire 
range 

St®  0.2  over  entire 
range 

3.8  x  105  <  Re<  1.14  x  106 

Near  wake  region  shrinks  considerably,  drag 
decreases  sharply,  periodic  vortex  shedding  is 
no  longer  ableto  bedetected  experimentally. 

Cd  *  0.08  at  the  critical 
point,  then  slowly 
rises  to  »  0.12 

Considerable 
uncertainty  exists 
in  the  literature 

3.  The  Fluent  Code 

Fluent  is  a  computational  fluid  dynamics  computer  code  developed  and  marketed  by  Fluent 
Inc.  [11].  The  code  solves  the  equations  for  conservation  of  mass,  momentum,  energy  and 
other  rel  evant  fl  ui  d  vari  abl  es  usi  ng  a  cel  I -centred  fi  nite-vol  u  me  method.  Fi  rstthefl  ui  d  domai  n 
is  divided  into  a  large  number  of  discrete  control  volumes  (also  known  as  cells)  using  a  pre¬ 
processor  code  which  creates  a  computational  mesh  on  which  the  equations  can  be  solved. 
The  mesh i  ng  softw are  avai  I  abl  e  w  i  th  F I  u ent  i  s  cal  I  ed  G ambi  t.  Th i  s  softw are  al  I  ow s  the  u se  of 
several  types  of  computational  cells  including  triangular,  quadrilateral,  hexahedral, 
tetrahedral,  pyramidal,  prismatic  and  hybrid  meshes. 

Once  the  fluid  domain  has  been  meshed,  the  governing  equations  (in  integral  form)  are 
applied  to  each  discrete  control  volume  and  used  to  construct  a  set  of  non-linear  algebraic 
equationsforthediscretedependent  variables.  Fluent  then  offers  theuser  a  number  of  choices 
for  the  algorithm  used  to  solvethese  equations,  including  coupled  explicit,  coupled  implicit, 
and  segregated  solvers.  In  all  the  calculations  reported  hereonly  the  segregated  solver  has 
been  used.  In  this  approach  the  governing  equations  are  solved  sequentially.  Si  nee  these 
equations  are  non-linear  they  are  first  linearised  using  an  implicit  method.  This  produces  a 
scalar  system  of  equations  containing  only  one  equation  per  computational  cell  perdegreeof 
freedom.  A  point  implicit  (Gauss-Siedel)  linear  equation  solver  is  then  used  in  conjunction 
with  an  algebraic  multi  grid  (A  MG)  method  to  solve  the  resultant  scalar  system  of  equations 
for  the  dependent  vari  abl  e  i  n  each  cel  I .  Si  nee  the  equati  ons  are  non-l  i  near  several  i  terati  ons  of 
the  solution  loop  must  be  performed  before  a  converged  solution  is  obtained. 
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3.1  Pressure-Velocity  Coupling 

Using  this  approach,  an  equation  for  each  component  of  the  momentum  equation  and  then 
thecontinuity  equation  aresolved  sequentially.  Once  thethree  components  of  velocity  have 
been  calculated  for  each  cell  using  this  sequential  system  the  velocities  may  not  satisfy  the 
continuity  equation  and  so  a  "Poisson-type”  equation  for  a  pressure  correction  is  derived  from 
thecontinuity  equation  and  the  linearised  momentum  equations.  This  pressure  correction 
equati  on  i  s  then  sol  ved  to  obtai  n  the  necessary  correcti  ons  to  the  pressu re  and  vel  ocity  fields 
such  that  continuity  is  satisfied. 

A I  though  the  pressure  vari  abl  e  appears  i  n  each  of  the  component  momentum  equati  ons  each 
of  these  equations  is  solved  by  treating  the  relevant  component  of  velocity  as  the  unknown 
variableand  thepressurefield  istaken  to  be  that  from  theprevious  iteration.  In  this  sequential 
procedurethe  continuity  equation  isused  as  an  equati  on  for  the  pressu  re.  However  pressure 
does  not  appear  explicitly  in  thecontinuity  equation  for  incompressi  bleflows  (which  are  the 
only  flows  considered  in  this  report)  and  so  a  procedure  must  be  devised  to  introduce 
pressure  into  this  equation.  Fluent  provides  methods  based  on  the  SIMPLE  (Semi-Implicit 
Method  for  Pressu  re- Linked  Equations)  family  of  algorithms  to  do  this. 

The  basi  c  SI  M  P  L  E  al  gori  th  m  uses  a  rel  ati  onsh  i  p  betw  een  vel  oci  ty  and  p  ressu  re  correcti  ons  to 
enforce  mass  conservation  and  to  obtain  thepressurefi  eld. TheSIM  PLEC  algorithm  (SIM  PL  E- 
Consi  stent)  is  a  variation  of  the  SIMPLE  algorithm  which  uses  a  more  refined  expression  for 
the  vari  able  flux  through  each  of  the  cell  faces.  ThePISO  pressure-velocity  coupling  scheme 
( P  ressu  r  e- 1  mp  I  i  ci  t  w  i  th  Sp  I  i  tti  n  g  of  O  perators)  i  s  al  so  p  art  of  th  e  SI  M  P  L  E  f  ami  I  y  of  al  gor  i  th  ms 
and  is  based  on  a  higher  degree  of  approximation  for  therel  ati  onshipbetw  een  thecorrecti  ons 
f  or  p  ressu  re  and  vel  oci  ty .  The  P I  SO  al  gori  th  m  takes  more  C  PU  ti  me  per  sol  ver  i  terati  on,  but  i  t 
can  dramatically  decrease  the  number  of  iterations  required  for  convergence,  especially  for 
transient  problems.  The  PISO  algorithm  also  allows  Fluent  to  obtain  solutions  on  highly 
skewed  meshes  in  approximately  the  same  number  of  iterations  as  required  for  more 
orthogonal  meshes. 

3.2  Spatial  Discretisation 

Fluent  stores  discrete  values  of  the  variables  at  the  cell  centres,  however  values  of  the 
variables  are  required  at  the  cell  faces  for  the  convection  terms  in  the  equations  and  these 
must  be  interpolated  from  the  cel  I  centre  values.  In  Fluent  6.1  this  is  accomplished  using  an 
upwind  scheme  and  there  are  several  alternatives:  first-order  upwind,  second-order  upwind, 
power  law,  and  QUICK. 

I  nthefirst-order  upwind  scheme  quantities  at  cell  faces  aredetermined  by  assuming  that  the 
cell-centre  values  of  any  variable  represent  a  cell-average  value  and  hold  throughout  the 
entire  cell.  Hence  the  face  quantities  are  identical  to  the  cell  quantities.  The  Power-Law 
scheme  i  s  an  i  improvement  whi  ch  i  nterpol  ates  theface  val  ue  of  a  vari  abl  e  by  usi  ng  the  exact 
solution  to  a  one-dimensional  convection-diffusion  equation.  When  theflow  isdominated  by 
convecti  on  thi  s  i  mpl  i  es  that  theface  val  ue  of  the  vari  abl  e  i  s  effectively  equal  to  the  cel  I  val  ue 
in  the  upwind  direction.  If  the  flow  is  weak  and  diffusion  stronger  then  the  Power-Law 
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scheme  amounts  to  a  simple  linear  average  of  the  value  of  the  variable  at  the  current  cell 
location  andtheupstreamceH.TheSecond-OrderUpwindSchemeprovidestruesecond  order 
accuracy  by  performing  a  Taylor  series  expansion  of  the  cel  I -centred  solution  about  the  cel  I 
centroid.  TheQU  I CK  scheme  also  al  lows  cal  culati  on  of  a  higher-order  value  of  the  converted 
variable  at  the  cell  face  by  using  a  weighted  average  of  second-order  upwind  and  central 
interpolation. 

Fluent  6.2  contains  enhanced  discretisation  schemeswhich  havebeen  motivated  by  the  new 
L  E  S  an  d  D  E  S  tu  r  bu  I  en  ce  mod  el  s  avai  I  abl  e  i  n  th  i  s  versi  on  of  the  cod  e.  Becau  se  accu  rate  sp  ati  al 
discretisation  is  crucial  in  LES  version  6.2  of  Fluent  contains  a  new  second-order  central 
differencing  scheme  for  the  discretisation  of  convective  terms.  Even  higher  order  upwind- 
biased  schemes  such  as  second-order  upwind,  QUICK  and  third-order  MUSCL  schemes 
introduceconsiderableamounts  of  numerical  diffusi on. This  is  particularly  detrimental  in  LES 
cal  cu  I  ati  ons  because  nu  meri  cal  d  iff u si  on,  no  matter  how  smal  I ,  can  easi  ly  overw  hel  m  physi  cal 
diffusion.  Central  differencing  schemes  are  well  known  to  haveeither  very  low,  or  sometimes 
zero,  numerical  diffusion,  and  hence  are  ideal  for  LES  calculations.  This  unfortunately  can 
lead  to  another  problem  however,  in  that  the  lack  of  any  numerical  damping  can  allow 
unphysical  oscillations  in  the  solutions  to  develop.  To  overcome  this  possibility  a  bounded 
central  differencing  scheme  was  also  introduced  in  version  6.2  of  Fluent.  This  scheme 
essenti  al  I  y  detects  any  osci  1 1  ati  ons  i  n  the  sol  uti  on  fi  el  d  w  i  th  a  wavel  ength  I  ess  than  tw  i  ce  the 
local  grid  spacing.  Itthen  suppresses  these  oscillations  by  switching  to  upwind  schemes  of 
varyi  ng  orders,  dependi  ng  on  the  severity  of  the  osci  1 1  ati  ons,  but  only  i  n  the  regi  ons  of  space 
where  these  osci  1 1  ati  ons  occur.  Because  of  these  features,  Fluent  version  6.2wasused  for  all 
the  LES  simulations  described  in  section  6. 

4.  T urbulence  M  odels 

Tu  rbu  I  ent  f  I  ow  s  are  characteri  sed  by  vel  oci  ty  f  i  el  d  s  w  h  i  ch  f  I  uctu  ate  rapi  d  I  y  both  i  n  space  and 
time.  Si  nee  these  fluctuations  occur  over  several  orders  of  magnitude  it  is  computationally 
very  expensive  to  construct  a  grid  which  directly  simulates  both  the  small  scale  and  high 
frequency  fl  uctu  ati  ons  for  probl  ems  of  practi  cal  engi  neeri  ng  si  gnifi  cance.  Two  methods  can  be 
used  to  eliminate  the  need  to  resolve  these  small  scales  and  high  frequencies:  Reynolds 
Averaging  and  Filtering. 

4.1  Reynolds  Averaging 

In  the  Reynolds  Averaged  approach  all  flow  variables  aredivided  into  a  mean  component 
(which  can  vary  slowly  in  time)  and  a  rapidly  fluctuating  component  and  then  all  equations 
are  time  averaged  to  remove  the  rapidly  fluctuating  components.  In  the  Navier- Stokes 
equati  on  theti  me  averagi  ng  i  ntrod  uces  new  terms  w  hi  ch  i  nvol  ve  mean  val  ues  of  prod  ucts  of 
rapi  d  I  y  varyi  ng  quanti  ti  es.  These  new  terms  are  know  n  as  the  Rey  nol  ds  Stresses,  and  sol  uti  on 
of  the  equations  initially  involves  the  construction  of  suitable  models  to  represent  these 
Reynolds  Stresses.  One  approach  is  to  treat  the  time  averaged  terms  as  additional  viscous 
stresses  produced  by  the  turbulence  in  the  flow.  In  the  Boussinesq  approximation  the 
Reynolds  Stresses  are  assumed  to  have  a  form  identical  to  the  viscous  stresses  in  the 
momentum  equation,  apart  from  a  multiplicative  term  known  as  the  turbulent  viscosity,  |_iT. 
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Fluent  provides  several  turbulence  models  based  on  the  Boussinesq  approach:  theSpalart- 
Allmaras  model  [45],thek-s  model  [24, 30,53],  and  the  k-co  model  [53].  As  discussed  in  section 
2  however,  noneof  these  model  shave  proven  to  be  parti  cu  I  arly  effective  at  capturing  theti  me 
dependent  features  of  theflow  past  a  sphereat  moderately  large  Reynolds  numbers.  This  is  in 
contrast  to  the  fi  nd  ings  of  I  accari  no  eta/.  [20],  who  carried  out  URA  NS  (Unsteady  Reynolds 
Averaged  Navier  Stokes)  simulations  of  flow  around  square  cylinders  and  a  wall-mounted 
cube  and  found  good  agreement  of  the  time-dependent  features  of  the  flow  (such  as  the 
frequency  of  vortex  shedding)  with  avail  able  experimental  data.  Constantinescu  and  Squires 
[7]  however  have  noted  that  U  RA  N  S  si  mul  ati  ons  past  spheres  appear  to  be  practi  cal  I  y  steady- 
state  when  compared  with  the  results  obtained  from  URA  NS  si  mul  ati  ons  past  cylinders,  and 
they  attri butethesedifferences to  thefact that thecoherence of  thewakestructures  in  theflow 
past  spheres  is  less  significant  as  compared  to  theflow  past  cylinders. 

Given  the  results  of  Drikakis[10]  and  Constantinescu  and  Squires  [7]  showing  that  URA  NS 
cal  cul  ati  ons  are  unabl  eto  capture theful  I  ti  me-dependent  structureof  fl  ow  past  a  sphere,  and 
the  results  of  Tomboulides  et  al.  [51],  Constantinescu  and  Squires  [7]  and  Kim  [25],  which 
demonstratesuperi or  predictivecapabilities  using  both  LES  (Large Eddy  Simulation)  and  DES 
(Detached  Eddy  Simulation)  schemes,  no  attempts  were  made  to  simulate  turbulent  flow 
structures  around  a  sphere  using  URANS  calculations.  All  turbulent  flow  calculations 
reported  here  used  the  LES  scheme,  which  is  described  inmoredetail  in  the  next  section. 

4.2  Filtering 

The  alternati  veto  Reynolds  averaging  is  filtering.  The  idea  behind  this  approach  is  to  filter  the 
time-dependent  Navier- Stokes  equation  in  either  Fourier  (wave-number)  space  or 
configuration  (physical)  space.  This  filtering  process  effectively  filters  out  turbulent  eddies 
whose  scales  are  smaller  than  thefi  Iter  width,  which  is  usually  taken  to  be  the  mesh  size.  A 
simulation  using  this  approach  is  known  as  a  Large  Eddy  Simulation  (LES).  As  with  Reynolds 
averaging  however,  the  filtering  process  creates  additional  unknown  terms  which  must  be 
modelled  in  order  to  providedosureto  theset  of  equations.  These  terms  are  the  sub-grid  scale 
stresses  and  Fluent  6.2  provides  several  modelsfor  these  stresses.  Thesimplest  of  theseisthe 
model  originally  proposed  by  Smagorinsky  [44],  in  which  thesub-grid  seal e stresses (SGS)  are 
computed  using  an  isotropiceddy  viscosity  approach  which  sets theSGSequal  totheproduct 
of  an  eddy- viscosity  and  the  resolved  rate  of  strain  tensor.  The  eddy  viscosity  is  then 
calculated  from  an  algebraic  expression  involving  the  product  of  a  model  constant  Cs,  the 
mod  ul  us  of  the  rate  of  strai  n  tensor,  and  an  expressi  on  i  nvol  vi  ng  thefi  I  ter  w  i  dth.  The  probl  em 
with  this  approach  is  that  there  is  no  single  value  of  the  constant  Cs  which  is  universally 
applicable  to  a  wide  range  of  flows.  Another  problem  is  that  the  model  is  not  applicable  to 
transitional  flows  because  the  basic  expression  always  gives  a  finite  value  for  the  SGS 
viscosity  even  in  laminar  regions  of  theflow,  as  long  as  there  is  a  velocity  gradient.  An 
additional  problem  is  that  an  ad  hoc  damping  is  needed  in  the  near-wall  region. 

Fluent  6.2  offers  several  LES  models  which  overcome  these  failings  of  the  original 
Smagorinsky  model.  The  WALE  model  (Wall  Adapting  Local  Eddy  Viscosity)  contains  a  near¬ 
wall  modification  in  theexpression  for  theSGS  viscosity  which  adapts  to  the  local  near-wall 
flow  structure  and  accounts  for  wall  damping  effects  without  using  a  damping  function 
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expl icitly.  Fluent  6.2  also  contains  the  Dynamic  Smagorinsky  Model  (DSM)  devised  by 
Germanoet  a/.  [12]  and  Lilly  [33].  In  this  procedureCs  is  dynamically  computed  during  the 
simulation  using  the  information  provided  by  the  smaller  scales  of  the  resolved  fields.  Cs 
d  eter  mi  n  ed  i  n  th  i  s  w  ay  var  i  es  w  i  th  ti  me  an  d  sp  ace  and  th  i  s  al  I  ow  s  th  e  Smago  r  i  n  sky  mod  el  to 
cope  with  transitional  flows  and  to  include  near-wall  damping  effects  in  a  natural  manner. 

Another  SGS  LES  model  in  Fluent  6.2  is  the  Localised  Dynamic  Kinetic  Energy  Model 
(LDKEM)  of  Kim  and  Menon  [28,29].  In  this  model  theSGSeddy  viscosity  is  computed  from 
an  expression  involving  thefilterwidth,theSGSturbulent  kinetic  energy,  and  a  (dynamically 
computed)  constant  Cs.  Thi  s  model  i  s  more  i  nvol  ved  i  n  that  an  additi  onal  transport  equati  on 
needs  to  be  sol  ved  for  the  SG S  tu  rbu  I  ent  ki  neti  c  energy,  but  the  L  D  K  E  M  has  several  d esi  rabl  e 
attributes  which  the  DSM  lacks.  In  particular,  using  the  SGS  turbulent  kinetic  energy  to 
parameteri sethe  SGS  stress  renders  the  LDKEM  better  suited  to  non-equilibrium  flows.  A 
problem  with  any  LES  calculation  is  the  treatment  of  wall  boundaries.  In  Fluent  6.2,  if  the 
mesh  isfineenough  to  resol vethelaminar  sublayer  (as  isthecasein  thesimulations discussed 
i  n  Secti  on  6)  the  use  of  w al  I  f  u ncti  ons  i  s  avoi  d ed  by  cal  cu I  ati  ng  the  w al  I  shear  stress  from  the 
laminar  stress-strain  relationship. 

Another  approach  is  known  as  Detached  Eddy  Simulation  (DES).  This  was  first  proposed  by 
Spalart[45, 46]  in  an  attempt  to  combine  the  most  favourable  aspects  of  RAN  Sand  LES.  DES 
reduces  to  a  RA NS  calculation  near  solid  boundaries  and  a  LES  calculation  away  from  the 
wall.  Fluent  6.2  offers  a  RANS/  LES  hybrid  model  based  on  the  Spalart-A  1 1  maras  turbulence 
model  nearthewall  and  a  one-equation  SGS  turbulence  model  away  from  the  wall  which 
reduces  to  an  algebraic  turbulent  viscosity  model  for  the  SGS  turbulence  far  from  the  wall. 


5.  Simulation  Results  for  Laminar  Flow 

5.1  Steady  axisymmetric  regime:  20  <  Re  <  210 

Section  2  describes  the  flow  over  a  sphere  as  being  both  axisymmetric  and  steady  up  to  a 
Reynolds  number  of  approximately  210.  To  test  the  software  in  this  range  a  run  was  first 
performed  at  Re  =  100  on  an  axisymmetric  mesh  containing  198,694  quadrilateral  cells.  The 
sphere  was  located  at  the  origin  and  had  a  diameter  of  1.0  m.  Theflow  wasintheX  direction 
with  a  velocity  of  0.147  cm/  s  and  the  viscosity  and  density  were  those  of  ambient  air,  i.e. 
1.7894x  10"5kg/  (m-s)  and  1.225kg/  m3  respectively.  The  mesh  extended  from  10m  upstream 
of  the  sphere  to  20  m  downstream  and  the  radial  extent  of  the  mesh  was  12  m.  The  laminar 
cal  cu  I  ati  on  w  as  ru  n  w  i  th  the  segregated  sol  ver  i  n  axi  sy  mmetri  c  geometry  usi  ng  the  cel  I  based 
gradient  option.  The  Pressure- Velocity  coupling  used  the  SIMPLE  method  and  the 
momentum  equations  were  solved  using  second-order  upwind  discretisation.  The  run  had 
converged  after  approximately  1000  iterations  and  thecalculated  valueofthedrag  coefficient 
wasCd  =1.087. 


In  order  to  check  any  dependence  of  this  result  on  the  resolution  of  the  grid  a  region  of  the 
mesh  between  X  =-4  m  to  X  =4  m  and  Y  =0to  Y  =4m  was  refined  using  the  Fluent  grid 
adaption  feature.  The  method  works  by  isotropically  subdividing  each  cell,  so  that  a 
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quadrilateral  cell  issplitinto  four  equivalent  quadrilateral  cel  Is.  This  increased  the  number  of 
cel  I  s  on  the  mesh  from  198  694  cel  I  s  to  240  470  cel  I  s.  The  val  u  e  of  the  d  rag  coeff  i  ci  ent  from  th  i  s 
si  mu  I  ati  on ,  to  f  ou  r  si  g  n  i  f  i  cant  f  i  g  u  res,  w  as  agai  n  C  d  =  1. 087. 

To  check  any  dependence  of  this  result  on  the  location  of  the  boundaries  a  new  mesh  was 
constructed  which  extended  from  15m  in  the  upstream  direction  to  25  m  in  the  downstream 
direction  and  out  to  15  m  in  the  radial  di  recti  on.  Thegrid  spacing  was  keptto  approximately  8 
mm  along  the  sphere  surface  and  the  height  of  the  first  cell  was  approximately  4  mm,  with 
successive  cell  heights  increasing  with  a  ratio  of  1.02.  The  total  number  of  cells  in  the  mesh 
w  as  277  344.  Th  e  cal  cu  I  ated  val  u  e  of  th  e  d  rag  coeff  i  ci  ent  w  as  C  d  =  1. 086.  Thisvaluewasalso 
confirmed  to  be  independent  of  grid  resolution  by  again  adapting  the  grid  in  the  region 
betweenX  =-4mtoX  =4mand  Y  =0toY  =4m.Thecalculated  value  of  the  drag  coefficient 
was  again  found  to  be  Cd  =  1.086. 

The  value  for  the  drag  coefficient  calculated  using  Fluent  on  the  smaller  of  the  unrefined 
meshes  (Cd  =  1.087)  compares  favourably  with  the  results  found  by  many  other  authors. 
Tabata  and  Itakura  [47]  have  performed  a  precise  computation  of  the  drag  coefficient  of  a 
sphere  in  the  range  between  Re  =  10  and  Re  =  200  using  two  finite  element  schemes  for 
axi  symmetric  flow  problems.  Oneisan  extension  of  the  standard  mixed  finite  element  method 
to  axi  symmetric  problems  while  the  other  is  an  extension  of  the  stabilised  finite  element 
method.  The  combination  of  these  two  methods  results  in  a  consistent  flux  method  which 
al  I  ows  them  to  obtai  n  error  esti  mates  for  the  d  rag  coeffi  ci  ent.  Thei  r  method  for  Re = 100  gi  ves 
Cd  =  1.0900  ±  0.0003.  Clift,  et  al.  [5]  cite  a  value  of  Cd  =  1.087,  while  Le  Clair,  et  al.  [31] 
calculated  a  value  of  Cd  =  1096  at  this  Reynolds  number.  M  ittal  [37]  has  used  a  Fourier- 
C  hebyshev  spectral  col  I  ocati  on  method  for  si  mul  ati  ng  fl  ow  past  spheres  and  spheroi  ds  i  n  the 
Reynolds  number  range  between  Re  =  50  and  Re  =  500  and  finds  Cd  =  1.09  at  Re  =  100. 

In  order  to  further  validate  the  software  in  this  Reynolds  number  range  a  series  of 
computations  was  performed  on  the  smaller,  unrefined  mesh,  between  Re=20and  Re=200. 
Figure  1  shows  the  calculated  values  of  the  drag  coefficient  in  this  range,  as  well  as  the 
experimental  correlation  of  Clift,  et  al.  [5]  and  Putnam  [41],  the  calculated  values  of  Tabata  and 
ltakura[47]  and  Mittal  [37]  and  the  experimental  results  of  Roos  and  Willmarth  [42].  Thereis 
excellent  agreement  between  the  Fluent  simulations  and  the  experimental  and  simulation 
results  of  other  authors. 
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Figure  1:  Comparison  of  computed  drag  coefficient  with  the  experimental  correlation  of  Clift,  et  ai. 

[5],  thecalculated  values  of  Tabata  and  Itakura  [47]  andM  ittal  [37],  andtheexperimental 
resu Its  of  R  oos  and  Willmarth  [42 ] 


5.2  U nsteady  planar-symmetric  regime:  Re  =  300 

Because  of  the  I  oss  of  axial  symmetry  at  higher  Reynolds  numbers  a  three-dimensional  grid  is 
needed  for  all  simulation  runs  conducted  above  Re =210.  This  was  constructed  bycreatinga 
sphere  of  lm  diameter  at  the  origin  and  then  enclosing  the  sphere  in  a  larger  sized  cube.  A 
smaller  cube  was  then  inserted  inside  the  sphere  and  then  six  hexagons  were  created  by 
connecting  the  corner  points  of  the  inner  cubes  to  the  corner  points  of  the  outer  cube.  All 
redundant  geometrical  entities  were  then  deleted  and  only  the  sphere  faces  retained.  The 
sphere  surface  was  then  meshed  with  quadrilaterals  having  side  lengths  of  approximately 
5  mm  and  then  hexagonal  cells  were  created  with  the  first  cell  height  being  approximately 
15  mm  and  further  cell  heights  having  expansion  ratios  of  approximately  1.08.  Thethree- 
di  mensi  onal  geometry  was  then  meshed  usi  ng  theGambit  mappi  ng  opti  on.  Theboundari  es  of 
the  mesh  extended  fromX  =-8mtoX  =10m,  Y  =-8mtoY  =8mand  Z  =-8mtoZ  =8m. 
The  total  number  of  hexahedral  cells  was  965 120. 

The  three-dimensional  grid  was  first  tested  by  running  a  simulation  at  Re  =  100.  The 
calculated  drag  coefficient  had  a  value  of  1.0979  and  the  value  of  the  lift  coefficient  was 
1.4480  x  10"7,  indicating  that  the  solution  is  effectively  axisymmetric  even  though  the 
cal  cu  I  ati  on  was  carri  ed  out  on  a  three-di  mensi  onal  mesh.  The  gri  d  was  then  refi  ned  usi  ng  gri  d 
adaptation  in  the  region  X  =-2  mto  X  =  2  m,  Y  =-2mtoY  =  2m  and  Z  =-2  m  to  Z  =  2  m. 
The  converged  value  of  the  drag  coefficient  was  then  Cd  =1.0985,  i.e.  approximately  0.05% 
higher.  This  indicates  that  the  mesh  has  sufficient  grid  resolution  near  the  surface  of  the 
sphere.  Theval  ue  1.0985  is  approximately  1%  higher  than  theval  ueof  1.087  calculated  on  the 
axi  sy  mmetri  c  mesh  though  andindi  cates  that  there  may  still  be  some  si  i  ght  d  ependence  of  the 
value  on  the  location  of  the  boundaries. 
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A  time-dependent  run  on  the  three-dimensional  mesh  was  then  performed  at  a  Reynolds 
number  of  300.  The  simulation  was  again  run  in  laminar  modeusing  the  segregated  solver, 
the  SIMPLE  Pressure- Velocity  coupling,  second-order  upwind  discretisation  and  the  cell 
based  gradient  option.  The  run  was  started  from  the  time-independent  solution  found  at 
Re  =  100  and  thesi  mul  ati  on  was  al  I  owed  to  run  without  perturbati  on  unti  I  numeri  cal  round¬ 
off  errors  resulted  in  the  appearance  of  a  regular  periodic  solution. 

It  was  found  that  a  considerable  period  of  time  was  required  to  elapse  before  the  motion 
settled  into  a  regular  periodic  pattern.  Figure  2  shows  a  plot  of  the  drag  coefficient  as  a 
function  of  time  between  t  =480  seconds  and  t  =560  seconds.  The  mean  drag  coefficient  and 
its  rms ampl itudecal cul ated  from  Fi gure  1  areO.668 and  0.004 respectively.  J ohnson  and  Patel 
[23]  obtained  0.656 and  0.004,  Constantinescu,  et  al.  [8]  also  obtained  0.656 and  0.003,  while 
Tomboulides,  et  al.  [50]  calculated  0.671  and  0.003.  Giacobello  [13]  calculated  a  mean  drag 
coeffi  ci  ent  of  0.658,  w  hi  I  e  K  i  m,  et  al .  [27]  cal  cul  ated  a  val  ue  of  0.657.  Theti  me  step  used  for  thi  s 
run  was  0.05  seconds.  Reduction  of  this  value  to  0.025  seconds  produced  no  change  in  the 
val  ue  of  the  mean  drag  coefficient. 

It  should  be  noted  that  Johnson  and  Patel  and  Constantinescu,  et  al.  used  a  very  similar 
meshing  scheme  and  that  both  had  an  effective  blockage  ratio,  defined  as  the  sphere's  frontal 
areadivided  by  the  computational  grid'sfrontal  area,  of  about  0.11%.  The  mesh  for  thethree- 
di  mensi  onal  cal  culations  descri  bed  here  has  a  bl  ockage  rati  o  approxi  mately  three ti  mes  this 
value,  i.e.  0.31%.  The  si  mu  I  ati  on  run  on  the  three-dimensional  mesh  at  Re  =  100  resulted  in  a 
drag  coeffi  ci  ent  approxi  matel  y  1%  hi  gher  than  the  val  ue  cal  cul  ated  on  the  two-di  mensi  onal 
axi  symmetric  mesh  (which  had  a  blockage  rati  oof  approximately  0.1%).  If  we  assume  that  a 
si  mi  I  ar  red  u  cti  on  i  n  the  val  u  e  of  the  d  rag  coeffi  ci  ent  w  ou  I  d  resu  1 1  f rom  a  F I  uent  ru  n  on  a  mesh 
having  a  blockage  ratio  of  approximately  0.1%  instead  of  0.3%  then  the  calculated  valueof 
0.668would  reduce  to  0.661,  which  is  less  than  onepercentdifferentfromtheresultobtained 
byjohnson  and  Patel  and  Constantinescu  and  Squires.  The  bl  ockage  ratiofor  thecal  culations 
of  Tomboulides,  eta/.  [50]  was  1.23%,  which  may  explain  their  higher  valueof  Cd  =0.671. 

The  Strouhal  number  calculated  from  Figure  2  is  0.133.  Johnson  and  Patel  [23]  obtained 
St  =0.137,  Constantinescu  et  al.  [8]  found  St  =0.136,  whileTomboulides,  et  al.  [50]  obtained 
St  =  0. 136.  G  i  acobel  I  o  [  13]  cal  cul  ated  a  val  ue  of  0. 134  w  hi  I  e  K  i  m,  et  al .  [27]  al  so  fou  nd  a  val  ue  of 
0.134.  The  value  measured  by  Sakomoto  and  Flaniu  [43]  at  Re  =  300  is  in  the  range  of 
0.15-0.16.  The  time  step  used  by  Johnson  and  Patel  [23]  was  0.05  seconds  while 
Tomboulides,  etal.  [50]  used  a  much  smaller  time  step  of  0.005  seconds.  Whilst  reducing  the 
time  step  in  our  calculation  from  0.05  seconds  to  0.025  seconds  produced  no  change  in  the 
mean  valueof  the  drag  coefficient,  it  did  result  in  a  slight  increasein  the  Strouhal  number, 
from  0.133  to  0.136.  This  compares  well  with  the  values  calculated  by  other  authors  and  the 
experimental  results. 
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Figure  2:  D  rag  coefficient  versus  time  at  Re  =300 


The  periodicity  of  theflow  at  Re =300  isdueto  the  regular  shedding  of  vorti  city.  Asdescri  bed 
by  Magarvey  and  Bishop  [34],  the  wake  consists  of  a  series  of  interconnected  vortex  loops 
above  Re  =  290,  while  Sakamoto  and  Haniu  [43]  show  that  hairpin  shaped  vortices  begin  to 
shed  periodically  at  Re  =  300.  To  observe  the  shedding  of  vorticity  in  the  simulation  results 
presented  here  we  have  used  the  Q  method  of  vortex  identification.  This  is  one  of  several 
methods  which  have  been  proposed  and  discussed  in  the  literature  recently.  Jeong  and 
H  ussain  [22]  conducted  an  extensive  series  of  tests  on  several  vortex  identification  methodsto 
determine  the  most  appropriate  method  of  identifying  a  vortex  in  a  turbulent  flow.  They 
considered  thefol lowing  methods: 

1.  Local  pressure  minimum 

2.  Closed  or  spiralling  streamlines  and  pathlines 

3.  Iso-contours  of  vorti  city  magnitude(i.e.  |©|,  where®  isthevorticity,co  =  V  x  u,and  u  is 
theflow  velocity) 

4.  The  occurrence  of  complex  eigenvalues  of  the  velocity  gradient  tensor  Vu,  which  is 
equi  val  ent  to  the  statement  that  the  d  i  scri  mi  nant  of  the  characteri  sti  c  equati  on  for  the 
ei  genval  ues  i  s  posi  ti  ve 

5.  TheQ  method,  whereQ  is  the  positive  second  invariant  of  the  deformation  tensor, 
i  .e.  Q  =  ( Qij  Qy  -  Sy  Sj )/  2,  where  fly  =  (8u  J  dx\  -  duj  dx\ )/  2  and 

Sj  =(5Ui/  dXj  +8Uj/  dx\ )/  2.  A  vortex  exists  when  Q  >  0 

6.  The  ^2  method,  where  a  vortex  is  defined  as  a  connected  region  of  space  with  two 
negative  eigenvalues  of  the  tensor  S2  +fi2  (whereS  and  Q.  aredefined  above).  Since 
S2  +fi2  is  symmetric  and  has  three  real  eigenvalues  A-i,  and  X3  this  definition  is 
equivalent  to  the  statement  that  X2  <  0  if  the  eigenvalues  are  ordered  such  that 

A,i  >  ^2  >  A.3 
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Their  conclusion  was  that  only  the  X2  method  represented  the  topology  and  geometry  of 
vortex  cores  correctly  for  the  I  arge  vari  ety  of  fl  ows  consi  dered  i  n  thei  r  anal  yses,  although  they 
did  also  state  that,  in  most  cases,  the  Q  and  X2  methods  did  result  in  similar  vortex  cores. 
Lesieur,  et  al.  [32]  have  also  done  a  comparison  of  vortex  identification  methods  including 
local  pressure  mini  mum,  iso-contours  of  vorti  city  magnitude,  positiveQ  and  negative^and 
thei  r  results  also  show  great  si  mi  I  ariti  es  between  vortex  cores  i  dentified  usi  ng  the  positive  Q 
and  negative^  methods.  Giacobello  [13]  has  recently  simulated  flow  around  a  sphere  at 
Re  =  300  using  a  spherical  polar  coordinate  system  and  a  Fourier-Chebyshev  collocated 
method  as  part  of  a  broader  study  looking  at  the  wake  structure  of  transversely  rotating 
spheres.  He  elucidated  the  vortex  structure  in  the  wake  using  iso-surfaces  of  streamwise 
vorticity,  iso-surfaces  of  vorticity  magnitude,  theQ  method,  and  the  ^2  method,  amongst 
others.  He  also  concluded  that  the  X2  method  provided  the  best  representation  of  the  vortex 
structure  in  the  flow,  although  he  also  noted  that  theQ  method  and  the^2  criteria  produced 
very  si  mi  lar  vortex  structures.  Given  that  thetwo  methods  produce  very  similar  structures  for 
thefl  ow  arou  nd  a  sphere,  and  that  the  Q  method  i  s  far  easi  er  to  i  mpl  ement  i  n  FI  uent,  w  e  have 
decided  to  use  theQ  method  for  the  simulation  results  described  here. 

Figure  3  shows  contours  of  velocity  magnitude  plotted  on  iso-contours  of  Q  at  the  time 
t  =  675.35  seconds.  G i  acobel  I  o  [  13]  noted  that  an  i  important  property  of  a  vortex  i  dentifi  cati  on 
criteria  (which  is  amply  satisfied  by  the^2  method)  is  that  it  should  be  insensitive  to  the 
threshold  valueusedforthevisualisation.  Wehavealsofound  thisto  betrueinthecaseofthe 
Q  method  as  quite  I  arge  changes  in  the  iso  value  lead  to  identical  vortex  structures.  Figure  3 
clearly  shows  the  shear  layer  separating  from  the  sphere  surface  and  the  tubular  vortical 
structure  i  n  the  far  wake. 

Figure  4  shows  a  different  perspective  of  the  same  plot.  This  clearly  shows  the  planar 
symmetry  which  still  persists  in  thefl  ow  atthis  Reynolds  number,  as  noted  by  Sakamoto  and 
Haniu  [43]  experimentally,  and  as  demonstrated  by  Mittal  [36]  in  his  numerical  simulations. 
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F i gu re  3:  C  on  tou rs  of  vel oci ty  magn itude at  Re  =  300  plotted  on  iso- con  tou rs  of  thesecon d  i n  vari an  t 
of  the  deformation  tensor 
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Jun  07,  2006 

FLUENT  6.2  (3d,  segregated,  lam,  unsteady) 


Figure  4:  Con  tours  of  velocity  magnitude  at  Re  =  300  plotted  on  iso-  con  tours  of  the  secon  d  i  n  vari  an  t 
of  the  deformation  tensor 
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6.  Simulation  Results  for  Turbulent  Flow 

As  discussed  insertion  2.5,  URA  NS  simulations  are  unable  to  capture  the  fine  scale  time 
depen  dent  features  of  the  turbulent  flow  around  aspheresuch  as  thehigh  frequency  (Kelvin- 
H  el  mholtz  instability)  modefound  in  the  experimental  results  of  Sakomoto  and  Haniu  [43],  or 
the  presence  of  the  higher  frequency  band  seen  in  the  simulations  by  Constantinescu  and 
Squires  [7].  Fortunately  Fluent  6.2  contains  both  LESand  DES  turbulence  models  as  well  as 
considerably  enhanced  discretisation  schemes  which  are  suited  to  the  effective  utilisation  of 
these  advanced  models. 

As  described  in  Section  4.2,  Fluent  6.2  offers  an  LES  model  with  the  original  Smagorinsky 
expression  for  the  sub-grid  scale  stresses  or  three  advanced  LES  models;  WALE,  LDKEM,  or 
the  Dynamic  Smagorinsky  model.  The  DES  model  in  Fluent  6.2  uses  a  Spalart-Allmaras 
turbulence  model  near  the  wall  and  a  one-equation  SGS  turbulence  model  away  from  the 
wall. 

Kim  [25]  has  used  theLES  model  on  unstructured  meshes  to  simulate  flow  around  a  sphere 
for  Reynolds  numbers  of  1.0 x  104  and  1.14 x  106  using  both  the  dynamic  Smagorinsky  model 
and  the  LDKEM,  butfound  thattherewaslittledifferenceintheresults  predicted  by  thetwo 
models.  Constantinescu,  et  al.  [8]  have  also  simulated  the  flow  over  a  sphere  at  a  Reynolds 
number  of  1.0  x  104  using  both  LES  and  DES  models  and  found  that  both  methods  gave 
essentially  the  same  results.  In  this  section  we  use  the  Fluent  LES  model  with  the  Dynamic 
Smagorinsky  SGS  model  to  simulate  flows  at  Re  =  1.0  x  104  and  Re  =  1.0  x  106.  Because  the 
boundary  layer  is  considerably  thinner  at  these  higher  Reynolds  numbers  we  increased  the 
grid  resolution  in  this  region  by  constructing  a  new  grid  containing  2.6  million  hexahedral 
cel  I  s.  The  method  of  construed  on  and  outer  I  i  mi  ts  of  the  mesh  remai  n  the  same  as  those  of  the 
earlier  mesh  but  the  cell  density  has  been  increased  over  the  surface  of  the  grid  and  in  the 
near-wake  region. 

6.1  Reynolds  number  =  104 

Theti  me-dependent  run  was  started  from  an  i  niti al  mean  fl ow  cal  cul ated  from  a  steady-state 
RA  NS  calculation  using  the  k-s  model.  The  inflow  velocity  was  0.146  m/  sand  the  time  step 
was  0.137  seconds.  Thi  s  corresponds  to  about  250  ti  mesteps  per  osci  1 1  ati  on  based  on  a  known 
experimental  Strouhal  number  of  approximately  0.2  at  this  Reynolds  number.  The  flow 
structure  at  Re  =  1.0  x  104  is  remarkably  different  to  that  at  Re  =  300  as  can  be  seen  from 
Figured  which  shows  the  vorti  city  structure  of  the  flow  using  theQ  method  described  inthe 
previous  section. 

The  additional  complexity  in  the  flow  is  also  evident  in  the  plot  of  drag  coefficient  as  a 
function  of  time,  shown  in  Figure  6.  In  addition  to  the  main  frequency  due  to  the  vortex 
shedding,whichoccursata  Strou  halnumberofO.  191  i  n  F  i  gu  re  6,  th  e  ti  me  h  i  story  al  so  sh  ow  s 
an  add  i ti  onal  hi  gh  frequency  component  i  n  the  range  1.65  to  1.82  w  hi  ch  i  s  associ  ated  w  i th  the 
Kelvin-FI  el  mholtz  instability  in  the  detached  shear  layers,  as  noted  by  Constantinescu  and 
Squires  [7].  These  results  agree  well  with  those  of  the  LES  simulations  of  Constantinescu,  et  al. 
[8],  who  found  a  value  of  St  =0.195  for  the  fundamental  shedding  frequency,  as  well  as  a 
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second  high  frequency  band  occurring  between  St  =  1.30  and  St  =  1.85.  Kim  [25]  found 
St  =0.182  using  an  unstructured  hybrid  mesh  of  2.46  mi  1 1  ion  cel  Is  using  LES  with  a  Dynamic 
Smagori  nsky  SGS  model .  The  mean  drag  coeffi  ci  ent  cal  cul  ated  over  the  I  ength  of  the  run  was 
0.387.  Constantinescu,  et  al.  [8]  obtained  0.393  ±  0.014  while  Kim  [25]  calculated  a  value  of 
0.438.  The  experi  mental  value  (from  Achenbach  [2])  isCd  =0.40±  0.01. 
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Contours  of  Velocity  Magnitude  (m/s)  (Time=6.Q492e-*Q3)  Jun  07,  2006 

FLUENT  6.2  (3d,  segregated,  LES,  unsteady) 


Figure5:  Contours  of  velocity  magnitude  at  Re  =  1.0  x  104  plotted  on  iso-contours  of  the  second 
invariant  of  the  deformation  tensor 

Figure  7  shows  the  Y  component  of  the  lateral  force  coefficient  as  a  function  of  time.  This 
displays  a  fundamentally  different  time  structure  to  that  shown  by  the  drag  coefficient.  The 
Strouhal  number  estimated  from  Figure7  lies  in  the  range  St  =0.151-  0.158  and  thereis  no 
evidence  of  a  high-frequency  component  associated  with  the  shear-layer  instability. 
Constantinescu  and  Squires  [7]  found  two  dominant  frequencies  in  the  ranges  St  =0.06-  0.08 
and  St  =0. 15  -  0. 16  and  al  so  noted  the  absence  of  a  hi  gh-frequency  component  associ  ated  w  i  th 
the  shear  I  ayer  i  nstabi  I  ity .  They  stated  that  thi  s  was  expl  ai  ned  by  thefact  that  the  vortex  tubes 
shed  at  thehigher  Strouhal  number  arepositioned  approximately  in  planes  perpendicular  to 
the  free  stream  di  recti  on  and  theresulti  ng  force  on  the  sphere  is  therefore  ori  ented  paral  I  el  to 
this  direction.  They  also  noted  that  the  difference  in  frequency  between  the  lateral  force 
oscillations  and  the  oscillations  in  the  drag  coefficient  also  occur  for  cylinders. 
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Figure  6:  D  rag  coefficient  versus  time  at  Re  =  1.0  xlO4 


time  (seconds) 


Figurel:  Y  component  of  lateral  force  coefficient  versus  time  at  Re  =  1.0  xlO4 

Figure  8  shows  streamlines  over  the  surface  of  the  sphere  from  two  different  viewpoints. 
These  cl  early  show  thelocation  of  boundary  layer  separation  as  well  as  the  extent  of  variation 
in  the  azimuthal  direction.  From  these  figures  we  estimate  that  separation  occurs  at  a  polar 
angleof  88.0°  ±1°.  This  is  slightly  larger  than  the  value  of  85.0°  ±  1°  found  by  Constantinescu, 
et  al.  [8]  in  their  LES  simulations  although  it  is  closer  to  the  value  found  by  Kim  [25],  who 
found  a  separati  on  angl  e  i  n  the  range  86°  ~  87°.  We  note  that  the  y+  val  ue  over  the  surface  of 
the  sphere  near  the  separati  on  regi  on  i  s  i  n  the  range  y+  =  0. 2  to  y+  =  0.7. 
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a  b 

Figure  8:  a)  Surface  streamlines  viewed  along  theY  axis  atRe  =  104.  b)  Surface  streamlines  viewed 
along  theZ  axis  at  Re  =  104. 


6.2  Reynolds  number  =  106 

To  perform  si  mulations  at  a  Reynolds  number  of  106thegrid  used  for  the  previous  simulation 
at  Re  =  104  was  further  refined  in  a  region  close  to  the  sphere  surface  to  ensure  adequate 
resolution  of  thethinner  boundary  layer  at  this  higher  Reynolds  number.  The  mesh  contained 
3.3 million  hexahedral  cellsand  resulted  in  awall  y+valuewhich  varied  between  0.3and  0.9 
over  the  surface  of  the  sphere  duri  ng  the  course  of  the  si  mulation. 

The  simulation  was  again  started  from  an  initial  mean  flow  calculated  from  a  steady-state 
RA  NS  calculation  using  thek-s  model.  The  inflow  velocity  was  14.6  m/  s  and  the  time  step 
used  was  0.001  seconds.  At  this  Reynolds  number  the  flow  is  well  above  the  "drag  crisis" 
which  occurs  at  approximately  Re  =  3.8x  lCf  and  wewould  therefore  expect  the  sphere  to 
exhibit  a  much  lower  drag  coefficient  and  that  separation  would  occur  much  closer  to  the 
downstream  stagnation  point.  This  is  qualitatively  evident  from  Figure  9,  which  shows  the 
vorticity  structure  in  theflow. 
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a 
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Figure9:  a)  Contoursof  velocity  magnitude  at  Re  =  1.0  xlO6  plotted  on  iso-contours  of  the  second 
invariant  of  the  deformation  tensor,  b)  Contours  of  velocity  magnitudeat  Re  =  1.0  xlO6 
plotted  on  iso-contours  of  the  second  invariant  of  thedeformation  tensor.  Enhanced  view  of 
Figure  9a. 


The  dramati  c  change  i  n  thefl ow  structure  i  s  easi  I y  seen  by  compari  ng  Fi  gure  5  and  Fi gure  9. 
The  close-up  view  in  Figure  9b  in  particular  shows  that  the  wake  at  the  higher  Reynolds 
number  is  much  narrower,  as  expected,  due  to  the  delayed  onset  of  flow  separation  in  the 
now  turbulent  boundary  layer. 

Figure  10  shows  the  drag  coeffi  dent  as  a  function  of  ti  me.  The  mean  val  ue  i  s  0.104.  Kim  [25] 
found  Cd  =  0.139 using  the dynamicSmagorinsky  model  andCd  =0.142  using  theLDKEM  at 
Re = 1. 14 x  106.  Constanti  nescu  and  Squ i  res  [ 7]  fou nd  Cd  =  0.080 at  Re = 1. 14 x  106  i  n  thei  r  D  ES 
simulation,  and  also  found  thatCd  =0.102  if  they  forced  the  boundary  layer  to  remain  fully 
turbulent  before  separation.  Wang  and  Kannan  [52]  used  an  overset  adaptive  cartesi  on/  prism 
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grid  method  with  2.1  million  cells  and  a  Spalart-Almaras  turbulence  model  and  found 
Cd  =  0.09.  J  i  ndal  ,etal.  [21]  al  so  conducted  an  L  ES  si  mul  ati  on  on  a  sphere  at  Re  =  1.14  x  106  on 
an  unstructured  tetrahedral  mesh  with  2.5  million  cells  using  the  dynamic  Smagorinsky 
model  and  the  instantaneous  logarithmic  law  ofthewall  to  model  thewall  shear  stress.  They 
found  Cd  =0.141.  The  experimental  valuefrom  Achenbach  [2]  isCd  =0.12. 


Figure  10:  D  rag  coefficient  versus  time  at  Re  =  1.0  xlO6 

Fi  gure  11  shows  streaml  i  nes  over  the  su  rface  of  the  sphere.  A I  though  there  i  s  obvi  ousl  y  more 
vari  ati  on  i  n  the  I  ocati  on  of  the  separati  on  I  i  ne  i  n  the  azi  muthal  d  i  recti  on  than  occurred  for  the 
case  of  Re  =1.0  x  104  i  t  i  s  sti  1 1  possi  bl  eto  esti  mate  an  average  separati  on  angl  e  of  121°  ±  2°.  Thi  s 
compares  well  with  the  results  of  Constantinescu  and  Squires  [7],  who  found  0s=12O°±2°at 
Re  =  1.14  x  106,  and  also  with  the  experimental  valuefrom  Achenbach  [2]  of  0S  =120°.  Wang 
and  Kannan  [52]  found  0S  =  115°,  Jindal,  et  al.  [21]  found  0S  ~  120°,  while  Kim  [25]  found 
0s~  100°.  As  noted  by  Kim  [25],  the  discrepancy  in  this  case  was  undoubtedly  due  to 
insufficient  resolution  in  the  boundary  layer. 

Figure  12  shows  the  frequency  spectrum  of  the  time  history  of  the  drag  coefficient.  Peak 
power  is  observed  at  St  =0.048  and  this  agrees  well  with  the  result  found  by  J  i  ndal ,  et  a/.  [21] 
who  found  a  peak  at  a  Strouhal  number  of  St~  0.05.  Constantinescu  and  Squires  [7]  however 
found  a  peak  i  n  the  frequency  spectrum  at  St  =  1.30  for  Re  =  1.14  x  106. 
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Figure  11:  a)  Surface  streamlines  viewed  aiong  theY  axis  atRe  =  106.  b)  Surface  streamlines  viewed 
along  theZ  axis  atRe  =  106. 


Figure  12:  Plot  of  frequency  spectrum  as  a  function  ofStrouhal  number  calculated  by  taking  theFFT 
of  Cd  versus  time  at  Re  =  1.0  xlO6 
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7.  Discussion  and  Conclusion 

In  this  section  we  briefly  review  the  simulation  results  for  the  various  flow  regimes.  A 
summary  of  the  data  for  the  computed  mean  drag  coefficient  (Cd),  Strouhal  number  (St)  and 
mean  separation  angle(0s)  isgiven  inTable2.  Between  Re=20and  Re=200,  wheretheflow  is 
steady-state,  axi symmetric  and  laminar  the  simulated  valuefor  the  drag  coefficient  was  in 
excellent  agreement  with  results  found  both  experimentally  and  computationally  by  other 
authors.  I  n  particular,  at  Re  =  100,  thedrag  coefficient  was  calculated  to  beCd  =1.086,  which 
compares  favourably  with  the  experimentally  calibrated  value  of  Cd  =  1087  found  by 
Clift, eta/.  [5]  and  the  precise  computation  of  Tabata  and  Itakura  [47],  which  found 
Cd  =1.0900  ±  0.0003. 

At  Re  =  300  the  flow  becomes  unsteady  and  transitions  into  periodic  time-dependent 
behavi  our  dueto  the  regular  shedding  of  hairpin  shaped  vortices.  The  mean  drag  coefficient, 
after  correction  for  blockage  ratio  effects,  was  calculated  to  have  the  value  0.661.  This 
compares  very  favourably  with  thevalueof  Cd  =0.656  obtained  by  both  Johnson  and  Patel 
[23]  and  Constantinescu,  etal.  [8],  Giacobello  [13]  calculated  Cd  =0.658,  while  Kim,  et  al.  [27] 
calculated  a  value  of  0.657.  The  corresponding  Strouhal  number  calculated  from  this 
simulation  was  St  =  0.133,  which  also  compares  favourably  with  the  value  found  by 
Giacobello  [13]  and  Kim,  et  al.  [27],  who  both  calculated  St  =0.134,  whilejohnson  and  Patel 
[23]  obtained  St  =0.137  and  Constantinescu,  et  al.  [8]  found  St  =0.136. 

At  Re  =  104 the  flow  in  the  wake  becomes  turbulent  while  the  flow  in  the  boundary  layer 
remains  laminar.  An  LES  simulation  using  the  Dynamic  Smagorinsky  sub-grid  scale  model 
with  a  fine  mesh  to  resolve  the  laminar  sublayer  was  used  for  this  Reynolds  number.  The 
simulated  results  were  in  good  agreementwith  both  experimental  results  and  other  numerical 
simulations.  The  time-dependent  behaviour  becomes  considerably  more  complex  and  in 
addition  to  the  mai  n  frequency  dueto  the  vortex  shedding,  which  now  occurs  at  a  Strouhal 
number  of  0.191,  the  time  hi  story  also  shows  an  additional  high  frequency  component  in  the 
range  1.65  to  1.82.  These  results  agree  well  with  those  of  the  LES  simulations  of 
Constantinescu,  et  al.  [8],  who  found  a  value  of  St  =  0.195  for  the  fundamental  shedding 
frequency,  as  well  asa  second  high  frequency  band  occurring  between  St  =1.30  and  St  =  1.85. 
The  mean  drag  coefficient  was  0.387,  which  agrees  with  the  value  obtained  by 
Constantinescu,  et  al.  [8]  of  0.393  ±  0.014  and  is  very  close  to  the  experimental  value  of 
Cd  =0.40±  0.01  (Achenbach  [2]).Streamline  plots  over  the  surface  of  the  sphere  showed  that 
boundary  I  ayer  separati  on  occurs  at  a  pol  ar  angl  e of  88.0°  ±1°,  w hi  ch  i  s  si  i  ghtl y  I  arger  than  the 
val  ue  of  85.0°  ±  l°found  by  Constantinescu,  et  al.  [8]  but  is  closer  to  the  value  calculated  by 
Kim  [25],  who  found  0sto  be  approximately  86°  ~  87°. 

At  a  Reynolds  number  of  approxi  mately  Re =3.8xicf  the  boundary  I  ayer  switches  from  being 
predomi  nantl  y  I  ami  nar  to  predomi  nantl  y  tu  rbul  ent,  hence  the  separati  on  poi  nt  moves  further 
downstream,  the  wake  narrows,  and  thedrag  coefficient  drops  considerably,  from  around 
Cd  =  0.40  at  Re  =  1.0  xlO4  to  approxi  mately  Cd  =0.08.  A  further  increase  in  Reynolds  number 
then  leadstoaslow  increase  in  thedrag  coefficient  again  until  Cd~0.20atRe~4.0xl06.The 
Fluent  LES  simulation  at  Re  =  1.0  xlO6  adequately  captured  this  behaviour.  The  vorticity 
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structure  of  the  flow  showed  that  the  wake  was  much  narrower  at  the  higher  Reynolds 
number.  The  average  separation  angle  was  found  to  be  121°  ±  2°,  which  agrees  with  the 
results  of  Constanti  nescu  and  Squi  res  [7],  w ho  found  0S  =  120°±  2°  at  Re = 1. 14  x  106,  and  al  so 
with  the experi mental  valuefrom  Achenbach  [2]  of  0S  =120°. 

Table2:  Summary  of  simulation  results  as  a  function  of  Reynolds  number 


Re 

cd 

(Fluent) 

Cd  (calc/expt) 

st 

(Fluent) 

st 

(calc/expt) 

6s 

(Fluent) 

0s 

(calc/expt) 

100 

1.087 

1.087  [5] 
1.0900  [47] 

1.09  [37] 

1.096  [31] 

Not  applicable 

Not  applicable 

Not 

applicable 

Not  applicable 

300 

0.661 

0.656  [23] 
0.656  [8] 

0.657  [27] 
0.658  [13] 
0.671  [50] 
0.671  [49] 

0.133 

0.134  [13] 

0.134  [27] 

0.136  [8] 

0.136  [50] 

0.137  [23] 

Not 

applicable 

Not  applicable 

1.0  xlO4 

0.387 

0.393  [8] 

0.438  [25] 
0.4Qt0.01  [2] 

0.191  (Stl) 

1.65-1.82  (St2) 

0.195  [8]  Stl 
0.181  [25]  Stl 
1.30-1.85[8] 

88.0°  ±1° 

85.0°  ±1°  [8] 
86.5°±1°  [25] 

1.0  xlO6 

0.104 

0.080  [7] 

0.09  [52] 

0.12  [2] 

0.139  [24] 
0.141  [21] 
0.142  [24] 

0.048 

0.05  [21] 

1.30  [7] 

121.0°  ±2° 

100°  [25] 

115° [52] 
120.0°  ±2°  [7] 
120°  [2] 

120°  [21] 

The  simulation  results  described  i  n  this  report  havei  1 1  ustratedthecapabi  I  ity  of  Fluent  version 
6. 2  to  accurately  reproduce  typical  flow  structures  observed  on  this  generic  bluff  body  flow 
for  both  time  independent  time-dependent  and  laminar/  turbulent  flow  regimes. 
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